# 标准大气查询工具V1.4
# 大气参数计算以美国标准大气US1976为基准
# 可：基于海拔高度，计算对应高度的大气参数
# 可：基于海拔高度和真空速，计算对应马赫数、表速、总温、总压
# 可：基于海拔高度和表速，计算对应真空速、马赫数、总温、总压
# 可：基于海拔高度和马赫数，计算对应真空速、表速、总温、总压
# 可：基于海拔高度和总温，计算对应真空速、马赫数、表速、总压
# 注意事项：
# 1. 86km以上高度的数据缺乏准确度，118km以上的大气数据无法计算
# 2. 基于总温计算时，总温不应低于给定高度下的大气静温
# 3. 表速计算暂未考虑可压缩性修正，所以非亚声速条件结果可能有误，有时间会改
# 问题反馈：B站、微博@方向舵舵机改
# https://gitee.com/rudpt/atmos-tool
import math

class atmosphere():
    def __init__(self):
        self.p   = 0.0
        self.t   = 0.0
        self.rho = 0.0
        self.a   = 0.0       #声速
        self.pa   = 101325.0 #海平面压强
        self.ta   = 288.15   #海平面静温
        self.rhoa = 1.225    #海平面密度
        self.gamma = 1.4     #比热比

def USSA1976(alt):
    atmos = atmosphere()
    Rearth  = 6369.0 #地球半径km
    GMR     = 34.163195 
    #数据标签
    htab = [0.0     , 11.0          , 20.0          , 32.0          , 47.0          , 51.0          , 71.0          , 84.852    ]
    ttab = [288.15  , 216.65        , 216.65        , 228.65        , 270.65        , 270.65        , 214.65        , 186.946   ]
    ptab = [1.0     , 2.233611e-1   , 5.403295e-2   , 8.5666784e-3  , 1.0945601e-3  , 6.6063531e-4  , 3.9046834e-5  , 3.68501e-6]
    gtab = [-6.5    , 0.0           , 1.0           , 2.8           , 0.0           , -2.8          , -2.0          , 0.0       ]
    h = alt*Rearth/(alt+Rearth) #转换为位势高度
    #搜索高度所在标签范围
    i = 0
    j = 7
    iter = 0
    while True:
        k = int((i+j)/2)
        if h<htab[k]:
            j=k
        else:
            i=k
        if (j-i==1):
            break
        iter += 1
        if (iter>=4):
            print("所述高度无解")
            atmos.p = -1**0.5
            return atmos
    #列中间变量
    tgrad = gtab[i]
    tbase = ttab[i]
    deltah = h-htab[i]
    tlocal = tbase+tgrad*deltah
    if tgrad == 0.0:                                       #静压倍率
        delta = ptab[i]*math.exp(-GMR*deltah/tbase)
    else:
        delta = ptab[i]*pow(tbase/tlocal,GMR/tgrad)
    theta = tlocal/ttab[0]                                 #静温倍率
    sigma = delta/theta                                    #密度倍率
    #返回结果
    atmos.rho = atmos.rhoa*sigma
    atmos.p = atmos.pa*delta
    atmos.t = atmos.ta*theta
    atmos.a = pow(atmos.p/atmos.rho*atmos.gamma,0.5)
    return atmos


#检查atmos结果是否有问题
def atmos_error(atmos):
    if isinstance(atmos.p,complex) or isinstance(atmos.t,complex) or isinstance(atmos.rho,complex):
        print("无法给出结果,是否输入了高于118km的高度?")
        return True
    else:
        return False

#总温计算    
def total_temperature(ma,t,gamma):
    t0 = t*(1+(gamma-1)/2*ma**2)
    return t0

#总压计算
def total_pressure(ma,p,gamma):
    p0 = p*(1+(gamma-1)/2*ma**2)**(gamma/(gamma-1))
    return p0

#真空速转表速
def IAS_from_TAS(rho0,rho,vel):
    vela = vel*(rho/rho0)**0.5
    return vela

#动态数据输出
def move_print(mode,ma,t0,p0,vel,vela): 
    if mode=='0':
        return
    line_print()
    print("马赫数Ma     : ",format_number(ma))
    print("真空速(km/h) : ",format_number(vel*3.6))
    print("真空速(m/s)  : ",format_number(vel))
    print("总温(K)      : ",format_number(t0))
    print("总压(Pa)     : ",format_number(p0))
    print("表速(km/h)   : ",format_number(vela*3.6))

#静态数据输出
def alt_print(mode,atmos):
    if mode=='0':
        line_print()
    else:
        print("--------↓仅由高度决定↓--------")
    print("声速(km/h)   : ",format_number(atmos.a*3.6))
    print("声速(m/s)    : ",format_number(atmos.a))
    print("静温(K)      : ",format_number(atmos.t))
    print("静压(Pa)     : ",format_number(atmos.p))
    print("密度(kg/m^3) : ",format_number(atmos.rho))
    line_print()

#格式化输出
def format_number(num, float_precision=3):
    if abs(num) > 1e7 or (abs(num) < 1e-2 and abs(num) > 1e-5):  #这里可以根据需求调整阈值
        #使用科学计数法
        formatted = f"{num:.{float_precision}e}"                 
    else:
        #不使用科学计数法
        formatted = f"{num:.{float_precision}f}".rstrip('0').rstrip('.')
    return formatted

def line_print():
    print("==============================")

if __name__ == '__main__':
    while True:
        #功能选择循环
        while True:
            print('===标准大气参数计算程序V1.4===')
            print('------参考标准: USSA1976------')
            try:
                mode = (input("选择计算模式: \n  0-仅高度\n  1-高度+真空速 \n  2-高度+表速 \n  3-高度+马赫数 \n  4-高度+总温\n输入c退出\n"))
                break
            except ValueError:
                print("输入无效，请重新输入")
        if mode=='0':
            #参数输入循环
            while True:
                #参数输入次循环
                while True:
                    alt = (input("输入高度（单位: km）:\n输入c返回上一级\n"))
                    try:
                        alt = float(alt)
                        break
                    except ValueError:
                        if alt=='c':
                            break
                        print("输入无效，请重新输入")
                if alt=='c':
                    line_print()
                    break
                atmos = USSA1976(alt)
                if atmos_error(atmos):
                    continue
                alt_print(mode,atmos)
        elif mode=='1':
            while True:
                while True:
                    try:
                        alt,vel = (input("输入高度和真空速，空格分隔（单位: km, km/h）\n两输入均为c返回上一级:\n").split())
                    except ValueError:
                        print("输入无效，请重新输入")
                        continue
                    try:
                        alt = float(alt)
                        vel = float(vel)*1000.0/3600.0
                        break
                    except ValueError:
                        if alt=='c':
                            break
                        print("输入无效，请重新输入")
                if alt=='c':
                    line_print()
                    break
                atmos = USSA1976(alt)
                if atmos_error(atmos):
                    continue
                vela = IAS_from_TAS(atmos.rhoa,atmos.rho,vel)
                ma = vel/atmos.a
                t0 = total_temperature(ma,atmos.t,atmos.gamma)
                p0 = total_pressure(ma,atmos.p,atmos.gamma)
                move_print(mode,ma,t0,p0,vel,vela)
                alt_print(mode,atmos)
        elif mode=='2':
            while True:
                while True:
                    try:
                        alt,vela = (input("输入高度和表速，空格分隔（单位: km, km/h）\n两输入均为c返回上一级:\n").split())
                    except ValueError:
                        print("输入无效，请重新输入")
                        continue
                    try:
                        alt = float(alt)
                        vela = float(vela)*1000.0/3600.0
                        break
                    except ValueError:
                        if alt=='c':
                            break
                        print("输入无效，请重新输入")
                if alt=='c':
                    line_print()
                    break
                atmos = USSA1976(alt)
                if atmos_error(atmos):
                    continue
                vela = IAS_from_TAS(atmos.rhoa,atmos.rho,vel)
                ma = vel/atmos.a
                t0 = total_temperature(ma,atmos.t,atmos.gamma)
                p0 = total_pressure(ma,atmos.p,atmos.gamma)
                move_print(mode,ma,t0,p0,vel,vela)
                alt_print(mode,atmos)
        elif mode=='3':
            while True:
                while True:
                    try:
                        alt,ma = (input("输入高度和马赫数，空格分隔（高度单位: km）:\n两输入均为c返回上一级:\n").split())
                    except ValueError:
                        print("输入无效，请重新输入")
                        continue
                    try:
                        alt = float(alt)
                        ma  = float(ma)
                        break
                    except:
                        if alt=='c':
                            break
                        print("输入无效，请重新输入")
                if alt=='c':
                    line_print()
                    break
                atmos = USSA1976(alt)
                if atmos_error(atmos):
                    continue
                vel = atmos.a*ma
                vela = IAS_from_TAS(atmos.rhoa,atmos.rho,vel)
                t0 = total_temperature(ma,atmos.t,atmos.gamma)
                p0 = total_pressure(ma,atmos.p,atmos.gamma)
                move_print(mode,ma,t0,p0,vel,vela)
                alt_print(mode,atmos)
        elif mode=='4':
            while True:
                while True:
                    try:
                        alt,t0 = (input("输入高度和总温，空格分隔（单位: km, K）:\n两输入均为c返回上一级:\n").split())
                    except ValueError:
                        print("输入无效，请重新输入")
                        continue
                    try:
                        alt = float(alt)
                        t0  = float(t0)
                        break
                    except:
                        if alt=='c':
                            break
                        print("输入无效，请重新输入")
                if alt=='c':
                    line_print()
                    break
                atmos = USSA1976(alt)
                if atmos_error(atmos):
                    continue
                if t0<atmos.t:
                    print("输入总温低于给定高度静温，请重新输入")
                    continue
                ma = t0/atmos.t-1
                ma = ma*2/(atmos.gamma-1)
                ma = ma**0.5
                vel = atmos.a*ma
                vela = IAS_from_TAS(atmos.rhoa,atmos.rho,vel)
                t0 = total_temperature(ma,atmos.t,atmos.gamma)
                p0 = total_pressure(ma,atmos.p,atmos.gamma)
                move_print(mode,ma,t0,p0,vel,vela)
                alt_print(mode,atmos)
        elif mode=='c':
            line_print()
            break
    

